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1. Introduction 


Neuroimaging is the branch of medicine whose purpose is to provide visual information 
about the structure and the anatomy of the brain. The main techniques in clinics are: Computed 
Tomography (CT), Magnetic Resonance Imaging (MRI), Diffuse Optical Tomography (DOT), Positron 
Emission Tomography (PET) and Single Photon Emission Computed Tomography (SPECT). 


CT scanning uses X-rays crossing the sample to image sections of the specimen in study. 
Specimen can be a living being, a part of if (e.g. the abdomen, a knee, the head, ...) or whatever 
non-living object. X-rays travel ballistically inside most of the materials (living tissue 
included), so measuring absorption of X-rays we can guess the composition of the sample 
we are measuring. Changing the direction of injection of the X-rays and merging absorption 
data coming from multiple directions makes a planar reconstruction of the examined section 
of the sample possible. CT is invasive in the sense that irradiates the patient with ionizing 
radiation. A little dose is given to the patient during a single CT session, anyway. CT scanning 
of the head is typically used to detect skull fractures, brain injuries, aneurysms, strokes, brain 
tumors and arteriovenous malformations in the brain. 


MRI is based on nuclear magnetic resonance principles. It uses a strong static magnetic field 
to align the nuclear magnetization of hydrogen atoms of water in the body and then radio 
frequency fields are generated to alter this magnetization alignment. Several coils mounted 
on the scanner are then able to detect the magnetic field produced by the altered hydrogen 
atom magnetization and to relate the recovery time of these short-lasting magnetic fields to 
the environment in which resonant hydrogen atoms lie. MRI, using non-ionizing radiation, 
is generally considered non-invasive for the patient and provides greater contrast between 
different soft tissues than CT does. Magnetic resonance images of the head are mostly used to 
detect brain tumors. 
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PET and SPECT are nuclear medicine techniques and they are used in cancer localization. A 
radionuclide is injected in the body via blood flow and then tracked with detectors placed 
around the patient. 


Functional neuroimaging is a particular way to perform medical neuroimaging focused 
more on functionality rather than just resolve anatomical features. The goal of 
functional neuroimaging is to detect spatial and temporal changes of activated areas of 
the brain, by measuring related physiological or physical features. The most common 
techniques are functional Magnetic Resonance Imaging (fMRI), Electroencephalography (EEG), 
Magnetoencephalography (MEG) and, recently, functional Near InfraRed Spectroscopy (fNIRS) and 
Diffuse Optical Tomography (DOT). PET and SPECT can be considered functional images as they 
use a radionuclide concentration to distinguish between healthy tissue and tumoral tissue. 


Usually functional neuroimaging techniques are divided into two main classes: those with 
high spatial resolution (~ 1 cm or less) and those with high time resolution (~ 1 s or less). 


e High spatial resolution and poor time resolution techniques (fMRI, PET and SPECT). 
e High time resolution and poor spatial resolution techniques (EEG and MEG). 


Having good time resolution and sufficient spatial resolution, DOT lays between the two 
classes and it is often considered the optimal compromise for functional neuroimaging 
studies. Spatial resolution in DOT is not an intrinsic feature, being the technique 
diffusion-limited. Nevertheless, in brain activation studies, integrating optical data and a 
priori anatomic data from MRI scans, brings to huge enhancements in spatial localization. 


2. Functional near infraRed spectroscopy (fNIRS) 


The simplest way to perform tissue oximetry using harmless electromagnetic radiation is to 
shine the tissue by means of a continuous beam of infrared light and to collect the re-emitted 
or transmitted light. 


Hemoglobin is the key: oxygenated and deoxygenated states have considerably different 
absorption spectra and it turns out that, in first approximation, oxyhemoglobin and 
deoxyhemoglobin also are the two main chromophores in tissues into the InfraRed window. 
Because the underlining hypothesis states the presence of two chromophores, we need to 
inject radiation at two different wavelengths, possibly chosen where the spectra have the 
greatest differences (see Fig. 1). This spectroscopic technique uses constant light intensity 
and is known in literature as Continuous Wave Near InfraRed Spectroscopy (CW-NIRS) The most 
diffused wavelengths for CW-NIRS are 690 nm and 820 nm. 


2.1 Absorption spectroscopy 


Spectroscopy is the study of the interaction between radiation and matter, in particular 
absorption spectroscopy measures the absorption of radiation, in a selected range of 
frequency, in radiation-matter interaction. It is largely employed in analytical chemistry to 
check for the presence of elements or substances, being the absorption spectrum a sort of 
fingerprint and so characteristic of the substances. The simplest application is the detection of 
the amount of the substance present. 
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Fig. 1. The absorption spectrum, that is the molar extinction coefficient vs. the wavelength, of 
the oxygenated (red) and deoxygenated (blue) states of the hemoglobin. I and Ip are the 
chosen wavelengths for Near InfraRed Spectroscopy measurements. 


When a radiation beam of known intensity Ip crosses a layer of known width L of a certain 
medium (chromophore), part of the radiation is absorbed by the medium and the rest is 
transmitted on the other side of the layer (Fig. 2). Modeling this phenomenon supposing 
that each infinitesimal layer absorbs an amount of radiation dI and that this absorption is 
proportional to radiation intensity, brings to the following mathematical model. Calling ji, the 
absorption coefficient and z the radiation beam direction, the radiation intensity I measured 
is modeled by the Lambert-Beer law: 


dI = —paldz (1) 

I dI L 
f T= f rdd 2) 
I = Ipe bel (3) 


The absorption coefficient is measured in m~! and has a straightforward interpretation: 
la = r is the mean free path a photon travels prior to being absorbed. Moreover, it can 
be decomposed into two distinct factors: the absorbent medium concentration and its molar 
extinction coefficient e: 


Ha = [cle (4) 
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Clear Medium 
(CW Light) 


Fig. 2. Absorption Spectroscopy in the most simple geometry. Radiation is injected in the 
sample from one side and radiation intensity is measure on the other side. The Lambert-Beer 
law then allows to compute the chromophore concentration. 


The extinction coefficient depends on the radiation wavelength: e = e(A) and so the 
absorption coefficient fig = Ha (A). Such a functional dependence constitutes the absorption 
spectrum of the substance considered’. When more than a chromophore is present, each 
contribute is summed: 


Ha = [c1] €1 + [co] €2 + [c3] €3 +--+ = 2 [ci] €; (5) 


Media for which this model holds are non-diffusing media, often called clear media and 
concentration measurements are easily conducted knowing the extinction coefficient of the 
chromophore and the length of the path traveled by the radiation. In the many chromophore 
experiment, performing multiple sessions, changing the radiation wavelength, allows to write 
a linear system and to obtain the unknown concentration of each chromophore. 


Performing the same experiment described before using different media, it is found that 
for a vast class of substances this model is not accurate, as it doesn’t take into account 
for the diffusion of radiation within the medium (Fig. 3(a)). Scattering events are due 
to discontinuities in dielectric properties of the medium and can be modeled as collisions 
between the photons of the radiation and particles in the medium. Many scattering models 
have been developed in physics: in classical electromagnetism by Rutherford, Thomson and 
Mie and in quantum physics by Compton and Raman [Feynman et al. (1964); Mie (1908); 
Rutherford (1911)]. A second important class of media then have to be described, media often 
called turbid media. Living tissues belong to this second class. 


The Lambert-Beer model is still valid for this media, but the absorption coefficient alone is no 
more enough to describe the intensity drop found. A second coefficient ps, called scattering 


1 In many cases the absorption spectrum is given with respect to the frequency v = -4, being c the speed 
of light and n the refraction index of the medium. 
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coefficient, is then introduced. The Lambert-Beer law it is written as follows: 
I = Ipe (Mat Hs)L (6) 


The scattering coefficient here defined has a similar interpretation as the absorption coefficient: 
l = 7 is the mean free path between two scattering events. Within this class of substances, 
absolute concentration measurements are possible only for thin samples (for thicker samples 
more complicated approaches are necessary, see section 2.3.1). Incidentally it can be noticed 
that this disturbing factor provides a helpful effect. Radiation scatters statistically in every 
direction, even backwards. It is then possible to collect radiation from the same side of the 
sample used for the radiation injection and to design experiments in reflectance geometry 
(Fig. 3(b)). 


Turbid Medium 
(CW Light) 


Turbid Medium 
(CW Light) 


Ío 


(a) Absorption Spectroscopy using a (b) When using turbid media experiments 
turbid medium does not allow to compute in reflection geometry are possible. 
the chromophore concentration within the Scattering provides a positive probability 
sample because of scattering contributes to have radiation traveling backwards 
to intensity attenuation. after encountering scattering centers in 


the medium. 


Fig. 3. Absorption spectrocopy in turbid medium in trasmission geometry (Fig. 3(a)) and 
reflection geometry (Fig. 3(b)). 


A different approach is to irradiate the sample with a non-constant light intensity, but it 
becomes necessary to develop a more general model for light traveling inside turbid media. 
This will be done in section 2.3. First a simpler way to separate scattering contributes is 
explained. 


In eq. 6 it has been seen that both absorption and scattering take part in the same way 
to radiation attenuation and that is not trivial to discriminate between the two distinct 
contributions. There’s at least a case in which a simple way the separation is easily done. 
Suppose to have a medium in which the chromophore concentration [c (t)] slightly varies 
in time and suppose to sample spectroscopic data in reflection geometry for a certain time 
interval. This case is not that far from reality as it may seem at first sight, for example in 
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biological tissues chromophores concentrations change in time because of blood flow. Then 
for each instant t, a Lambert-Beer equation can be written as follow: 


I (tn) = Ine” {etn leDE+G) ' 


where scattering contribution usL has been grouped into the term G and the differential 
pathlength factor D has been introduced to take into account that photon travel paths longer 
than the source-detector separation because they penetrate inside the medium, scatters 
multiple times and then a little part reaches the detector traveling backwards. The differential 
pathlength factor depends only on the mean refractive index of the crossed medium. 
Hypothesizing that G doesn’t vary in time and is constant for small variations of chromophore 
concentrations and introducing the physical quantity absorbance A = In i = —([c]eDL+G), 
we can compute the variation in chromophore concentration between measurement collected 
at time t, and time tm in the following way: 


AAmn = Am — An (8) 
= —cmEDL — G + cneDL + G (9) 
= — (cm —¢n)eDL (10) 
= —ACmnEDL (11) 


The scattering contributes vanishes and the variation in chromophore concentration depends 
only on the two intensity measurements collected: 


AAmn 
A = — 12 
Cmn CDL (12) 
EPL E 13 
~ DL oe ua 
. 1 In lin 
= DL (in T In = (14) 
1 K 
= ZDI In (=) (15) 


Most of times a reference measurement is conducted and all the concentration values are then 
computed as variation from the reference value. 


2.2 Continuous wave NIRS (CW-NIRS) 


Hemodynamic activity of the brain is related through metabolic processes to the electrical 
activity of the neurons, then, performing hemoglobin concentration measurements on the 
cortical tissues, it is possible to obtain information about spatial localization and temporal 
behavior of neuronal activity. The spectroscopic tools described in section 2.1 provide a 
straightforward way to do this. It is possible to inject radiation into the tissues and collect 
the photons coming out, for example using a reflectance geometry. The usefulness of 
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this technique would be to measure what happens within the blood flow, neglecting the 
contributes from the surrounding environment. 


Collecting data at different times and then computing the absorbance variation values AA), 
and AA}, it is possible to write the following system of equations: 


AA, = (Eup, A [Hb] + epo,,, A [HbO2]) Da, L a6) 
AAj, = (Egba A [Hb] + eno, A [HbO2]) Da, L 
In matrix form it becomes: 
MAn, 
Dat _ fee | A [Hb] | 17) 
Aan | ~ 
Dit EHb,Az €HbO>,A1} LA [HbO2| 
and finally: 
~ii AA), 
| A [Hb] | 7 i | De (18) 
= AA 
A [HbO | EHb,A> EHbO>,A, Dit 


Being this a two equations systems with the two unknowns A [Hb] and A [HbOy], a solution 
to the linear system is possible and the oxyhemoglobin and deoxyhemoglobin concentration 
variations with respect to a reference value can be computed (Eq. 18). If data are sampled 
for an interval of time, solving the linear system for each time sample, the time course of the 
concentrations is given. 


Derived quantities such as total hemoglobin [tHb] are often plotted vs. time in order to 
improve the data visualization. 


A [tHb] = A [Hb] + A [HbO,] (19) 


An important limitation of CW-NIRS is the possibility to obtain just concentration variations. 
Having absolute concentration values, hemoglobin saturation sO could be calculated: 
[HbO;] 


sO = T) (20) 


In principle there is no depth sensitivity with a single source-detector pair in CW-NIRS. 
Information on the composition of the media crossed by the CW light can be measured as 
an average. To increase light penetration large source-detector separation have been used (4-5 
cm) with the obvious consequence of increasing the sampling volume and thus worsening 
the spatial resolution. Increasing the complexity of the system, the use of multi-source and 
multi-detector tomographic or topographic arrangements could overcome this limitation and 
better depth sensitivity [Boas et al. (2004); Strangman et al. (2002)]. 


Laboratory prototypes and commercial instrumentations have been developed through the 
years to measure oxygenation of the muscles and of the cortical regions of the brain. A 
comprehensive list can be found in Wolf et al. (2007). 
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2.3 Photon migration approach 


A different way to deal with diffusion and discriminate between scattering and absorption 
contributes is to use non-constant light intensities in the injected radiation. A possible 
way to describe the physics of the problem is to start from the electric and magnetic fields 
associated with the radiation and to develop a rigorous mathematical theory, taking into 
account scattering, diffraction and interference. This is quite complicated, especially because 
of multiple scattering. 


A more practical physical model of photon propagation into diffuse media based on energy 
flow is the Radiative Transfer Equation (RTE). This model was first developed to describe 
the non-interacting neutrons diffusion in nuclear reactors [Sanchez & McCormick (1982)]. 
The hypothesis of non-interacting particles reasonably applies to photons in multiple elastic 
scattering regime. The RTE model neglects polarization effects. This is not a problem though, 
because even if light coherence exists, this is going to be completely lost after a few scattering 
events. 


2.3.1 The radiative transfer equation 


The average power that at position r and time t flows through the unit area oriented in the 
direction of the unit vector §, due to photons within a unit frequency band centered at v, that 
are moving within the unit solid angle around § is the spectral radiance I; (r,8,t,v). Thus at 
time t through the unit area dÈ lying in r, oriented as 8, within the unit solid angle dO and in 
the frequency interval [v,v + dv] flows a power dP given by: 


dP = I, (r,8,t,v) |p» 8|dZdQdv (21) 


Energy density ae is simply related to spectral radiance. The energy dE per unit frequency 
and per unit solid angle that crossed the area dÈ oriented as 8, in the time interval dt is: 


dE [,dEdt k 


dV dXvdt v e 


Where v is the speed of radiation in the medium, v = £. Energy density is the number of 
photons and so I; is proportional to the number of photons in the unit volume, with frequency 
v, that at time t are moving in the direction 8. As we deal with media in which the radiation 
frequency doesn’t change during propagation (elastic scattering), we can integrate I; over the 
frequency v and obtain the radiance I (r, &,t). 


Integrating also over the solid angle and dividing by v gives the quantity u(r,t) = 


I(r,8,t)dQ : = z ops . 
larai, measured in |Jm~°|, that represents the energy density at position r and time 


t. Finally the photon density can be obtained dividing the energy density by the energy of a 
single photon: 

u (r,t) 
hv 
The RTE is an integro-differential equation stating the balance of the incoming and outgoing 
radiation along the direction 8, at the time t, inside the volume element dV at the position r 


n(r,t) = (23) 
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[Ishimaru (1978); Martelli et al. (2010)]: 


g Te VI- (pat os) E+ Hs f p(8,8)I(r,8, t) dO! + 24) 
v ot Ar 


Each term appearing has a precise significance: 


e iol _, temporal change of energy. 


e —8. VI = change due to energy flow. 
e — (ua + ps) I > energy drop due to absorption and scattering. 


A al 


© tus fan p (8, 8’) I(r, 8’, t) dO! => energy gain from radiation coming from every direction, 
but scattered into the direction &, being p (8, 8’) the phase function, that is the probability 
that a photon flowing in the direction & is scattered into the direction 8’. 


e g=4q(r,8,t) = Radiation source inside dV [Wm Ssr~1] ; 
The RTE has two immediate properties: 


1. If Ip is a solution for a non-absorbing medium with time impulsive source term q (r, 8,t) = 
go (r, 8) 6 (t), then if the absorption coefficient is iq the solution is: 


I = Ige he? (25) 
This property is a generalization of the Lambert-Beer law [Sassaroli & Fantini (2004); 


Tsuchiya (2001)]. 


2. If I is the Green function for a medium having extinction coefficient wp = Ma + Ms and 
albedo a = ne then for a medium having extinction coefficient u¥ and the same albedo, 
the solution is scaled in this way: 


He (26) 


This scaling property is known as the similarity principle [Zege et al. (1991)]. 


Because of the high complexity of the RTE, no analytical solutions are available and 
approximation methods are then applied. Numerical methods like the Finite Difference Method 
or the Finite Elements Method (FEM) require a large amount of computational power, but are 
often applied [Arridge (1995); Arridge et al. (2000); Arridge & Hebden (1997); Arridge & 
Schweiger (1995); Arridge et al. (1993)]. Stochastic methods like the Monte Carlo are widely 
used in many biological applications [Boas et al. (2002); Fang & Boas (2009)]. 


2.3.2 The Px approximation 


Many ways to simplify the RTE have been developed through the years, but most methods 
are based on the Py approximation. Two derived quantities of interest are the photon fluence ® 
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and the photon flux J: 
® (r,t) = | I(r, 8,t) dQ (27) 
4r 
J(r,t) = f I (r, 8, t) êdQ (28) 
4r 


Radiance and source term are expanded in spherical harmonics [Boas (1996); Jackson (1999)]: 


+œ | 241 

1st E yE orn (r,t) Yim (8) (29) 
1=0 m=-—l1 
+oo | ] 

grs DE E fai (rt) Yin (8) (30) 
1=0 m=—1 


Phase function is expanded in spherical harmonics too, but first the reasonable hypothesis 
that the scattering amplitude is only dependent on the change in direction of the photon is 
made. With this assumption the probability of scattering from a direction § into a direction 8’ 
depends only on the angle between & and 4’: 


p(é-s)=)> iz iP; (8: 8’) (31) 


I 
M 
Mm 
>Q 

Ea 


(8°) Yim (8) (32) 


Where P; is the Legendre polynomial of order / and the angular addition rule has been used 
[Jackson (1999)]. The normalization factors 4/ arl and ael are introduced for convenience. 


Truncating the series at the ] = N term brings to the Py approximation of the RTE [Arridge 


(1999)]. In the P, approximation the radiance?, source and phase function reduce to: 
I(r,8,t) = G(r) + oe ir t)- 8 (33) 
4r 4r 
g(r, 8,t) = q0 (r,t) + Ža (r,t) -8 (34) 
p(s-8) = 80+ pees (35) 


Where we used the definition of fluence 27 and flux 28, obtaining: 


D (r,t) = poo (r,t) (36) 
; 5 (1-1 (r,t) — $11 (7,4) 
= (1-1 (r,t) + 1, (r,t)) (37) 
1,0 (r,t) 


J (r,t) = 


e3 


2 As will be realized soon, this brings to assume the radiance to be nearly isotropic. This is a 
good approximation for biological tissues and in general for high-albedo media, because of the few 
absorption events relative to the scattering events. 
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For the source term approximation the quantities monopole moment qo and dipole moment 
qi have been introduced: 


qo (r,t) = qoo (r,t) (38) 
5 (H1,- 1(r,t) — qi, (r,t)) 
qı (r,t) = wa (M1, 1 (r,t) +911 (7, #)) (39) 
gio (r,t) 


Normalizing gq to unity, g1 is the average cosine of the scattering angle: 


go = 1 (40) 
g= <8-s'> (41) 


Finally, substituting these approximated expressions into 24, we obtain: 


oP (r 
1P OD L uD (r,t) Y: J (r,t) + go (r,t) (a2) 
A ee 5V® (r,t) + a1 (r,t) wa 


Where we made use of the definition of reduced scattering coefficient y4: 


Us = (1 — g1) Ms (44) 


2.3.3 The diffusion approximation 


In a high-albedo (predominantly scattering) medium, the time for a substantial change in flux 
J is much longer than the time to traverse one transport mean free path. Thus, over one 
transport mean free path itself, the fractional change in flux is much less than unity. The 
diffusion approximation results from making the following assumptions: 


OJ (r,t) 
at =0 (45) 


qı (r,t) =0 (46) 


Dropping the dipole moment of the source is reasonable assuming an isotropic source. Eq. 43 
then becomes Fick's law: 


1 
J (r,t E Ne (r,t) (47) 
ee) 
Usually the factor multiplying the gradient is called diffusion coefficient, D = Gay We 


notice that in high-albedo medium (u4 >> pq) it doesn’t depend on absorption and can be 


calculated as D ~ mE 


3 This is consistent with the similarity principle previously introduced [Martelli et al. (2010)]. 
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Substituting into Eq. 42 we finally have the diffusion equation (DE): 


1 0® (r,t) 
v ot 


= —ua® (r,t) + V-DV® (r,t) + qo (r,t) (48) 


The diffusion equation models many problems in Physics and Chemistry, such as heat 
propagation, ion diffusion and semiconductor doping. Analytical solution to the DE in photon 
migration have been carried out in symmetrical geometries such as the infinite medium, the 
semi-infinite medium or the semi-infinite medium with a high-absorbent inclusion. 


In order to solve the DE, appropriate boundary conditions have to be declared. The most 
general boundary condition is the Partial Current Boundary Condition (PCBC): 


(r,t) +2ADV®(r,t)-aA=0 reav (49) 


Where A is a coefficient depending on the refractive index and is due to Fresnel reflection 
[Contini et al. (1997); Martelli et al. (2010)]. In practical applications, though, the most 
commonly applied boundary conditions are the Extrapolated Boundary Condition (EBC) and 
the Zero Boundary Condition (ZBC): 


r © OVext (50) 
reoVv (51) 


Where dVext is a surface distant 2AD from the boundary (extrapolated surface). 


CW, TR and FD solutions for simple geometries have been calculated [Martelli et al. (2010)]. 


2.3.4 FD-NIRS 


Frequency Domain Near InfraRed Spectroscopy (FD-NIRS) makes use of intensity modulated 
laser light (typically at radio frequencies), injecting it into the sample. The remitted wave 
is demodulated to obtain amplitude and phase as a function of the modulation frequency. 
The tissue acts like a low-pass filter: the amplitude is a decreasing function of the frequency 
and the phase typically increases. The analytical expressions for amplitude and phase can 
be obtained by a Fourier transform of the TR theoretical expression. Estimation of optical 
properties can thus be performed using the photon migration theory (see Fig. 4). 


Some commercial device using frequency domain techniques is available and measurements 
and basic research often used this approach to study biological tissues oxygenation [D’Arceuil 
et al. (2005); Fantini et al. (1995)]. 


2.3.5 TR-NIRS 


In section 2.2 a way to measure oxygenated and deoxygenated hemoglobin concentrations 
has been described. Using absorption spectroscopy techniques it is possible to investigate in 
vivo the oxygenation status of the superficial tissues and in particular of the cortical areas of 
the brain. Thus, CW-NIRS is a useful non-invasive technique for functional neuroimaging. 
As previously seen, the issue underlying Near InfraRed Spectroscopy is scattering. The 
way CW-NIRS deals with scattering is to perform a basal measurement and then to refer all 
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I (w) R (p,w) 
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Fig. 4. Scheme of reflection geometry FD-NIRS measurement. Intensity modulated light 
I (w) travels in the medium and part of it can be collected by a detector at a source-detector 
distance p. Amplitude and phase of the collected light R (p, w) depend on pig and pf. 


I (t) R (p,t) 


\ [N— 


Fig. 5. Scheme of reflection geometry TD-NIRS measurement. Injected light travels in the 
medium and part of it can be collected by a detector at a source-detector distance p. Knowing 
the temporal shape of the injected light I (t), the temporal shape of the collected light R (p, t) 
depends on pi, and y4. 


successive data to the baseline values, making the assumption that scattering is constant in 
time. This makes concentration variations the only possible obtainable data. 


In section 2.3.1 the physical model for radiation propagation inside tissues has been 
introduced. It will be now shown how RTE gives a practically useful tool to perform a kind of 
measurement in which absolute hemoglobin concentrations values can be sampled. 


In many biological applications the semi-infinite geometry is assumed and light intensity 
measurement are conducted in reflectance geometry (see Fig. 5). In this case the Green 
function of the Diffusion Equation with the Extrapolated Boundary Condition, expressed as 
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the detected power per unit area (usually known as reflectance R), takes the form: 


1 ie pot z 2 
575 e Dot e Ha | zoe Do — Zpe Do (52) 
2 (4r Dv)” 4 £5/2 


R (p,t) = 


Where p is the source-detector distance, zo = mi Zp = Zo + 2ze and ze = 2AD. 


Looking at the the Green function of the diffusion equation it is possible to notice that 
absorption and scattering contributes are naturally separated in the model. Plotting these 
solutions for different values of the absorption coefficient and of the scattering coefficient 
evidence this property (see Fig. 6). This technique based on the injection of light with a known 
temporal shape (typical a pulse-like shape) and the estimation of the optical properties from 
the emitted light temporal shape is called Time-Resolved Near InfraRed Spectroscopy (TR-NIRS). 


The light is collected by means of photomultipliers tubes (PMT) or avalanche photodiodes 
(APD) usually in Time Correlated Single Photon Counting (TCSPC) regime. Detailed 
descriptions of the way light is collected by light detectors is given in [Donati (1998)] and 
of TCSPC is given in [O'Connor & Phillips (1984)] and [Becker (2005)]. 


Estimation of the optical properties can be performed from all the analytical expressions 
for the TR response of a diffusive medium. Absorption and scattering coefficients can 
be computed by means of an inversion algorithm [Press et al. (1992)]. The instrumental 
response function (IRF) due to temporal dispersion in fiber, temporal jitter of detectors and 
finite pulse width of light sources has to be taken into account. Due to the linearity and 
time-invariance of the transport problem, the detected response is the convolution between 
the IRF and the theoretical response of the medium. Thus, a Levenberg-Marquardt non-linear 
minimization between experimental data and theoretical curve convoluted with the IRF is 
generally performed. 


An estimation of the optical properties of a diffusive medium can be also performed by means 
of Monte Carlo simulations. A certain number of simulations are needed, then by using the 
scaling properties of the RTE it is possible to implement a fast search of the correct simulation 
[Pifferi et al. (1998)]. These two method enable an absolute estimation of optical properties. 


A simple estimation of the variation of the absorption coefficient can be easily performed if 
only absorption changes are assumed (and usually are). If we collect at different times two 
reflectance curves Ro (p, t) and Rj (p,t), it is straightforward that: 


—_ 1 Rı (p,#) 
M= an (Rey) m 


In principle TR techniques provides a richer insight than CW to the problem of non-invasively 
probing of a diffusive medium. These approaches can discriminate between absorption 
and scattering contributions and derive absolute values for the hemodynamic parameters 
[Patterson et al. (1989)]. This however can be obtained only in simple geometries like the 
infinite or the semi-infinite homogeneous models. In a real heterogeneous medium like the 
human tissues and in particular the human head it is easier to derive changes with respect to a 
baseline or effective average parameters rather than absolute values. Advanced time-resolved 
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perturbation models for more complicated geometries, like multi-layered media, have been 
derived in the last years, but they would require the use of a priori information (e.g., the 
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(b) Dependance from absorption (jig = 0.04, 0.08, 0.16, 0.32 cm7!, u; = 20 
cm~!). 


Fig. 6. Normalized reflectance from the semi-infinite model, with source-detector distance 
p=2 cm. Reflectance vs. time is plotted in semi-logarithmic scale. 
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anatomy of the head as provided by a MRI scan) for their practical and effective use [Martelli 
et al. (2005)]. 


The actual potentiality of time-resolved techniques relies on an easier approach to the problem 
of depth sensitivity. As mentioned before, to enhance depth sensitivity, CW systems use 
large source-detector distance and multi-distance approaches. Conversely, depth sensitivity 
in TR-NIRS can be intrinsic in the measured reflectance curve having information on the time 
of flight of the photons in tissues. 


Near InfraRed Spectroscopy applications to the brain spread from rehabilitation monitoring 
to neural plasticity studies, from localization of cortical activation in motor or cognitive tasks 
to stroke diagnosis, and more. Motor tasks such as finger tapping or hand grasping are easy 
to perform and hundreds of studies have been published [Torricelli et al. (2007)]. Cognitive 
tasks present more difficulties in activation interpretation, nevertheless a number of study is 
available in literature [Bandettini et al. (1997); Butti et al. (2009); Heekeren et al. (1997)]. Fig. 7 
and Fig. 8 show two examples of TR-NIRS data. 


Fig. 7. A [Hb] (blue) and A [HbO3] (red) concentrations (in arbitrary units) vs. time (in 
minutes) collected from the frontal area of the brain during a cognitive task of sustained 
attention. Vertical lines indicated the beginning and the end of the task. The average 
concentration value for the first two minutes block was used as baseline and then subtracted 
to all data. 


3. Diffuse Optical Tomography (DOT) 


Since the very beginning of the history of NIRS, many efforts have been made to improve 
space resolution and data accuracy. Near InfraRed Spectroscopy provides functional 
information about the oxygenation status of the explored tissues, but a single source-detector 
pair is only able to probe a small underlying area. In this section a way to perform functional 
neuroimaging starting form NIRS data will be introduced. Production of two dimensional 
maps, often called topography, by linear interpolation of NIRS data was the first application to 
be investigated. This imaging technique is usually called Diffuse Optical Imaging (DOI) and, 
adding depth sensitivity, is known as Diffuse Optical Tomography (DOT). 
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Fig. 8. Oxyhemoglobin (in red) and deoxyhemoglobin (in blue) concentrations for a hand 
grasping task. Single-subject, right hand movement, average of 10 repetitions, TR-NIRS 
device, contrast enhanced for the deep layers (see section 3.1). 


To properly investigate a large area (larger than the common source-detector distance, that 
is 1-4 cm) a multi-channel approach is needed. The starting point is the arrangement of a 
number of sources and detectors to cover the area of interest and the management of the 
source-detector pairs. Depending on the physical dimensions of the area even a huge number 
could be arranged. No theoretical limitations on the number of optodes exists, however, 
actual technology features instrumentations with up to 32 sources and 32 detectors (CW) 
and 16 sources and 16 detection channels (TR). Shining NIR light on the skin and collecting 
the back-scattered light in a reflectance geometry from multiple points easily allows to draw 
a map of the hemoglobin concentrations around the area of interest. Spatial resolution is 
poor (depending on source-detector distances) and a little depth sensitivity (associated to the 
overlapping measurements in CW or intrinsic in TR) is obtained. 


Concentrations time series can be plotted vs time or a spatial map can be built (see Fig. 9). 
A movie-like map can be obtained if each frame corresponds to a time sample. More often, 
especially in books and publications, averaged concentrations data over time intervals are 
used to obtain maps relative to a time block of seconds or minutes, just as the example in Fig 
9. 
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Fig. 9. The concentration time series can be visualized as spatial map. Each time sample 
belongs to an average position between the associated source-detector pair. The map is built 
by linear interpolation with the neighboring points concentration values. 


More complicated ways to image optical data have been developed. Using MRI anatomical 
information and using suitable values for the average optical properties of the tissues a inverse 
problem can be established. Localized alterations of the optical properties, corresponding to a 
increased absorption and related to the hemodynamic responses of neural activations, can be 
put in correspondence with cortical features and a better localization of brain activations can 
be performed. 


DOT systems can consist of little more than a probe with fiber-optic sources and detectors, 
a piece of dedicated hardware about the size of a small suitcase and a laptop computer. 
Systems can be much larger, depending primarily on the type of laser source and detectors 
employed, but the approach generally offers a degree of portability unobtainable with many 
other modalities. For this reason, looking into the future, DOT may be ideally suited for 
clinical applications such as bedside monitoring of cerebral oxygenation. 


3.1 Time resolved DOT 


A fundamental point in NIRS measurements, whatever technique we are using, is the ability to 
separate systemic hemodynamic changes occurring in the superficial tissues, such as the skin, 
from functional hemodynamic changes related to brain activations. In order to reach this goal 
is fundamental to discuss NIRS depth sensitivity. Depth sensitivity is not intrinsic in CW- and 
TR-NIRS measurements, but can be achieved using proper optodes configurations or using 
time of flight information. Multi-distance source-detector approaches, both with CW [Saager 
& Berger (2005)] and with TR [Liebert et al. (2004)] techniques, have been proposed to improve 
depth selectivity and sensitivity. Single-distance approaches have also been discussed [Selb 
et al. (2005); Steinbrink et al. (2001)]. Contini et al. (2007) proposed a different approach to add 
depth information to the tissue probing problem, based on time-domain contrast functions. 
Wabnitz et al. (2008) discussed a method for depth selectivity based on time windows and 
moments of time-of-flight distributions for TD-NIRS. 
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Time-resolved curves statistically contain information about tissues which light photons pass 
through. Photons reaching the detector at early times surely have been back-scattered from the 
superficial tissues and thus can bring information about the superficial layers only. Photons 
reaching the detector at late times have travelled inside the tissues for a longer time and 
statistically could have probed deeper layers and carry information about them. We usually 
refer to them as early photons and late photons. 


A simple approach to the depth probing problem could be to develop a contrast function, 
considering just two domains: superficial layers and deep layers. 


In order to develop models to discriminate between the variations of the absorption coefficient 
in the superficial layers (Av?) and of the variation of the absorption coefficient in deep 
layers (AuPOWN), the quantity time-dependent photon path length, usually called Mean 
time-dependent Path Length (MPL) is introduced [Steinbrink et al. (2001)]. As its own name 
suggests, MPL is the average length of the path travelled by photons in a specific layer and 
can be calculated from the mean time of flight of the photons. 


The reflectance curve is divided into time intervals Tt, called time gates (TG), and the total 
counts (usually referred as intensity) are calculated for each time gatet. 


Being Ip (t) the intensity for the non-absorbent medium, LU? (t) and LPOWN (t) the mean 
time-dependent path lengths in the superficial layers and in the deep layers, respectively, the 
expression for the intensity I (t) in a general absorbent medium can be expressed as: 


I (t) = i (t) go (Aug? LYP ()+AppOWNLPOWN (4) ) (54) 


In the simplest case, it is possible to think that early photons carry information only about the 
superficial layers and that late photons are not affected by a superficial inhomogeneity, and 
carry information only about the deep layers. Thus, we can write for the absorption coefficient 
changes: 


UP __ 1 I (Te) 
Mi = -T 9 (Frey) i 
DOWN _ 1 I(t) 
Aha IPOWN (qij) in (;, =) an 


where Te is the mean time of flight of photons in a early time gate and 7 is the mean time 
of flight for photons in a late time gate. Despite giving a general idea about how to obtain 
information about contributes of different layers, this model is pretty rough, indeed. The 
hypothesis of late photons not affected by superficial layers has to be rejected and more 
complicated expressions are often necessary to perform depth selection [Contini (2007)]. It 
can be proved that late photons carry both information about superficial and deep layers and 
that superficial layer absorption variations affect the deep layers absorption estimation. The 
subtraction of the superficial contribute is thus necessary in the estimation of AuPOWN., Such 


te: 
1 The total counts are the time integral of the reflectance curve: Igate (t) = fi?” R (t) dt 
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a subtraction can be performed as follow: 


UP 1 I (Te) 
^ha = ~TUP Cz) (ie) Pa 


DOWN _ 1 Iq) I(t) 
Ha = —IDOWN (n) ( The). & a Pa 


Oxy- and deoxy-hemoglobin concentrations for the upper layer and the lower layer can be 
easily obtained via Lambert-Beer law, from the estimated A4. An example of application of 
this stratigraphic model is shown in Fig. 10, where data where collected during a Valsalva 
maneuver, that forces an increase of the systemic blood volume, resulting in an increase of 
both oxyhemoglobin and deoxyhemoglobin concentrations in the superficial layers of the 
skin. 


e — AHbO, DOWN 


Fig. 10. Estimated AHbO, and AHb in the superficial layer (UP, continuous line) and in the 
brain (DOWN, dotted line) during the Valsalva maneuver without normalization. The two 
vertical dashed lines indicate the beginning and the end of the task period, respectively. 
During the Valsalva maneuver the sistemic blood volume dramatically increases, while the 
local blood volume in the brain is less affected by the maneuver effects. A proper depth 
selectivity remarks this differences. 


3.2 Multimodality approach 


Recently a great interest in multimodality approaches has grown. Merging the advantages of 
different imaging and functional techniques, such as co-registration of NIRS with MRI [Merritt 
et al. (2002)], blood flow monitors, fMRI [Torricelli et al. (2007)] and PET, gives the possibility 
to build anatomo-functional images and movies to largely improve information visualization. 
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Moreover, the comparison with standard clinical techniques such EEG or MEG can lead to a 
clinical standardization of the NIRS signal and push Near InfraRed Spectroscopy towards a 
regular clinical use [nEUROPt Project (2008-2012)]. 


Using anatomical magnetic resonance priors to perform MRI-guided optical reconstruction 
dramatically improves the spatial resolution of the diffuse optical tomography techniques 
[Boas & Dale (2005)]. Simulations of light propagation are run into the MRI head volume, 
using Monte Carlo statistical methods such as Fang & Boas (2009), and then the cortical 
activation profile is obtained solving an inverse problem. For a more detailed description 
of it see Caffini et al. (2011). 


In the event of an unavailability of the subject-specific anatomy the efficient use of an 
MRI anatomical atlas has been demonstrated [Caffini et al. (2010)]. Fig. 11 reports an 
atlas reconstruction of a cortical activation during a visual protocol of pattern reversal 
checkerboard. 


Fig. 11. Atlas map visualization of an MRI-guided optical reconstruction of the brain 
activation in the visual cortex (the brain is seen from the back), during a pattern reversal 
checkerboard protocol. 


4. Conclusions 


Near InfraRed Spectroscopy applied in vivo to cortical tissues has been widely investigated 
through the last two decades and big steps have been done. 


From the technical point of view, the laser technology, and in particular the introduction of 
semiconductor laser diodes, has helped to fabricate compact and clinical instrumentations. 
Moreover, the improvements in light detectors has permitted to develop accurate devices, for 
example, in the time-resolved technology, the red-extended photocathodes well increased the 
quantum efficiency in the near red region. 


Nevertheless, a lot more has to be done. Pulsed lasers sources are getting better and better, 
especially concerning time stability and output power. Photomultiplier tubes is an efficient 
and well established technology, but, in NIRS application, suffers the need of high voltage 


www.intechopen.com 


72 Advances in Brain Imaging 


supplies and the necessity to work in a dark environment, to avoid background light. Single 
Photon Avalanche Diodes (SPAD) are the available technology that best fits the needs of Near 
InfraRed Spectroscopy. Actual work is to integrate SPADs into NIRS setups by means of 
increased light sensitive area and better quantum efficiency in the near red spectrum. For 
more information about SPAD detectors see [Cova et al. (2010)]. 


An interesting future perspective, is the so-called null-distance measurement setup, that is the 
possibility to collect light from the very same point of injection. Only time resolved techniques 
allow this possibility and a few successful efforts have been made in this direction [Pifferi et al. 
(2008)]. 


From the medical point of view, hundreds of physiological studies and psychological tasks 
have been carried out and a massive literature is available on the subject. For these reason, a 
larger clinical use of non-invasive optical imaging in the next years is expected. In particular, 
we expect time resolved NIRS to be the best candidate for this purpose. The four-year 
nEUROPt Project [nEUROPt Project (2008-2012)], financed by the European Union under The 
Seventh Framework Programme for research and technological development (FP7) for the 
period 2008-2011, and coordinated by the Authors, aims at the development and clinical 
validation of advanced non-invasive optical methodologies for in-vivo diagnosis, monitoring, 
and prognosis of major neurological diseases (stroke, epilepsy, ischemia), based on diffuse 
optical imaging by pulsed near infrared light. The consortium plans major developments in 
technology and data analysis that will enhance TD-NIRS with respect to spatial resolution, 
sensitivity, robustness of quantification as well as performance of related instruments in 
clinical diagnosis and monitoring. A strong clinical basis is being produced and the diagnostic 
value of TD-NIRS applications to brain study will be assessed, by putting using standard 
methodologies (such EEG) and new optical methods side by side, in a co-registration setup. 
The potential commercialization of TD-NIRS systems will be then evaluated by European 
system manufacturers. 
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